

*****
***** BEWARE - CODE IS COMPUTATIONALLY INTENSIVE
***** This code will run around 24hrs on a fast laptop/desktop computer.
***** If you simply want to assess replicability, reduce the number of bootstrap 
***** samples with the option reps(*), which currently stand at 5000.
***** Note, however, that small numbers will produce somewhat less precise CIs
*****


*estimate model-based survival estimates with bootstrap confidence intervals

set seed 52349085
use book-basedata-replication, clear

*estimate model to retreave sample for calculations of in-sample mean
stcox manip2 nomanip2 prevcris2 viol2 crisdur2 jointdem victory2 contig2 ///
	if _t<3650, strata(order5) cluster(dyadno) efron nohr texp(_t) tvc(manip2 prevcris2)

*calculate mean of each variable
foreach var of varlist med2 prevcris2 viol2 crisdur2 jointdem victory2 contig2 {
quietly sum `var' if e(sample)
local s_`var'=r(mean)
display `s_`var''
}

*predict survival with manipulative mediation at means
bsurvci if _t<3650, id(id) at(manip2 1 nomanip2 0 prevcris2 `s_prevcris2' ///
	viol2 `s_viol2' crisdur2 `s_crisdur2' jointdem `s_jointdem' ///
	victory2 `s_victory2' contig2 `s_contig2') strata(order5) ///
	ties(efron) texp(_t) tvc(manip2 prevcris2) gen(s_manip) replace reps(5000)

*predict survival with no manipulative mediation at means
bsurvci if _t<3650, id(id) at(manip2 0 nomanip2 1 prevcris2 `s_prevcris2' ///
	viol2 `s_viol2' crisdur2 `s_crisdur2' jointdem `s_jointdem' ///
	victory2 `s_victory2' contig2 `s_contig2') strata(order5) ///
	ties(efron) texp(_t) tvc(manip2 prevcris2) gen(s_nomanip) replace reps(5000)

*predict survival without mediation at means
bsurvci if _t<3650, id(id) at(manip2 0 nomanip2 0 prevcris2 `s_prevcris2' ///
	viol2 `s_viol2' crisdur2 `s_crisdur2' jointdem `s_jointdem' ///
	victory2 `s_victory2' contig2 `s_contig2') strata(order5) ///
	ties(efron) texp(_t) tvc(manip2 prevcris2) gen(s_none) replace reps(5000)

save bsurvci_results.dta, replace
	
*plot results for each stratum
local t=1
label var _t "Analysis time in days"
label var s_none`t' "No mediation"	
label var s_manip`t' "Manipulative mediation"	
drop s_nomanip*
twoway 	rarea s_none`t'_lb s_none`t'_ub _t if _t<3650, sort c(J J) lcol(gs10%30)  fcol(gs10%30) || ///
	rarea s_manip`t'_lb s_manip`t'_ub _t if _t<3650, sort c(J J) lcol(gs10%30)  fcol(gs10%30)  || /*rarea s_nomanip`t'_lb s_nomanip`t'_ub _t || */ ///
	line s_*`t' _t if _t<3650, c(J J J) lc(black black) lp(solid dash) ///
	legend(col(1) label(1 95%-CI) label(2 95%-CI) order(3 4) region(lwidth(none)) ) ///
	saving(strata`t', replace) t2title(`t' previous crises) sort 
graph export strata`t'.pdf, replace
forvalues t=2/5 {
	label var s_none`t' "No mediation"	
	label var s_manip`t' "Manipulative mediation"	
	*label var s_nomanip`t' "Non-manipulative mediation"	
	twoway 	rarea s_none`t'_lb s_none`t'_ub _t if _t<3650, sort c(J J) col(gs10%30)  fcol(gs10%30) || ///
		rarea s_manip`t'_lb s_manip`t'_ub _t if _t<3650, sort c(J J) col(gs10%30)  fcol(gs10%30)  || /*rarea s_nomanip`t'_lb s_nomanip`t'_ub _t || */ ///
		line s_*`t' _t if _t<3650, c(J J J) lc(black black) lp(solid dash) ///
		legend(col(1) label(1 95%-CI) label(2 95%-CI) order(3 4) region(lwidth(none))) ///
		saving(strata`t', replace) t2title(`t' previous crises) sort
	graph export strata`t'.pdf, replace
}
grc1leg strata1.gph strata2.gph strata3.gph strata4.gph strata5.gph, ///
	ycommon row(3) leg(strata1.gph) saving(full, replace) ///
	note(Note: Shaded areas depict 95% bootstrap CI.)
graph display, ysize(8)
graph export full.pdf, replace
graph export full.png, replace
